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During  1999-2000,  13  bottom  mounted  acoustic  Doppler  current  profilers  (ADCPs) 
and  12  wave/tide  gauges  were  deployed  along  two  lines  across  the  Korea/Tsushima 
Strait,  providing  long-term  measurements  of  currents  and  bottom  pressure.  Tidally 
analyzed  velocity  and  pressure  data  from  the  moorings  are  used  in  conjunction  with 
other  moored  ADCPs,  coastal  tide  gauge  measurements,  and  altimeter  measurements 
in  a  linear  barotropic  data  assimilation  model.  The  model  fits  the  vertically  averaged 
data  to  the  linear  shallow  water  equations  in  a  least-squares  sense  by  only  adjusting 
the  incoming  gravity  waves  along  the  boundaries.  Model  predictions  are  made  for  the 
Ol5  P15  Kl5  /i2,  N2,  M2,  S2,  and  K2  tides.  An  extensive  analysis  of  the  accuracy  of  the  M2 
surface-height  predictions  suggests  that  for  broad  regions  near  the  mooring  lines  and 
in  the  Jeju  Strait  the  amplitude  prediction  errors  are  less  than  0.5  cm.  Elsewhere,  the 
analysis  suggests  that  errors  range  from  1  to  4  cm  with  the  exception  of  small  regions 
where  the  tides  are  not  well  determined  by  the  dataset.  The  errors  in  the  model  pre¬ 
dictions  are  primarily  caused  by  bias  error  in  the  model’s  physics,  numerics,  and/or 
parameterization  as  opposed  to  random  errors  in  the  observational  data.  In  the  model 
predictions,  the  highest  ranges  in  sea  level  height  occur  for  tidal  constituents  M2,  S2, 
Kj,  Oj,  and  N2,  with  the  highest  magnitudes  of  tidal  velocities  occurring  for  M2,  Kl9 
S2,  and  Ov  The  tides  exhibit  a  complex  structure  in  which  diurnal  constituents  have 
higher  currents  relative  to  their  sea  level  height  ranges  than  semi-diurnal  constitu¬ 
ents. 
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1.  Introduction 

The  Korea/Tsushima  Strait  joins  the  East  China  Sea 
and  the  Yellow  Sea  (both  areas  of  strong  tides)  to  the  Ja¬ 
pan/East  Sea§  (an  area  of  weak  tides).  Ogura  (1933)  made 
the  first  comprehensive  study  of  the  tides  in  the  strait  using 
coastal  station  data  and  provided  a  basic  understanding 
of  the  tidal  heights,  tidal  dynamics,  and  amphidromic 
points.  Since  this  study,  many  others  have  investigated 
the  tides  in  the  area,  making  use  of  both  numerical  mod- 
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els  (Fang  and  Yang,  1988;  Choi  et  aL ,  1989;  Kang  et  al. , 
1991;  Kantha  et  al ,  1996)  and  observations  (Odamaki, 
1989;  Kang  et  al .,  1995).  These  studies  all  agree  in  a  quali¬ 
tative  sense  and  show  the  structure  of  the  tidal  heights  in 
the  strait.  Semidiurnal  tidal  height  amplitudes  peak  along 
the  south  coast  of  Korea,  with  a  general  trend  from  high 
values  in  the  East  China  Sea  to  low  values  in  the  Japan/ 
East  Sea;  amphidromic  points  for  these  constituents  oc¬ 
cur  just  north  of  the  strait  in  the  Japan/East  Sea  proper. 
Semidiurnal  tidal  height  phases  remain  constant  through 
much  of  the  strait,  but  increase  rapidly  from  west  to  east 
across  the  northern  boundary.  Diurnal  tidal  height  ampli¬ 
tudes  range  from  high  values  in  the  East  China  Sea  to 
low  values  in  the  Japan/East  Sea;  amphidromic  points  for 
these  constituents  occur  in  the  northwestern  comer  of  the 
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Fig.  1.  Locations  of  observational  moorings  (triangles),  coastal  tide  stations  (circles),  and  TOPEX/POSEIDON  altimeter  meas¬ 
urements  (diamonds).  Bathymetric  contours  for  the  50,  100,  200,  500,  and  1000  m  isobaths  are  calculated  from  the  SKKU 
dataset  and  drawn  as  thin  lines.  The  cities  of  Busan  (Pusan)  and  Fukuoka  (squares)  are  shown  for  reference.  Dark  shading 
indicates  land  areas.  Light  shading  indicates  ocean  areas  not  considered  by  the  data  assimilation  model.  Jeju  (Cheju)  Island  is 
the  large  island  to  the  south  of  Korea  in  the  far  western  portion  of  the  map.  Jeju  Strait  is  the  strait  between  Jeju  Island  and 
Korea.  Tsushima  Island  is  the  large  pair  of  islands  in  the  center  of  the  map.  The  Gotou  Islands  are  the  cluster  of  islands  in  the 
southern  portion  of  the  map  to  the  west  of  the  main  Japanese  island  of  Kyushu.  Honshu  is  the  main  Japanese  island  in  the 
northeast  portion  of  the  map. 


strait,  close  to  the  Korean  coast.  Diurnal  tidal  height 
phases  increase  along  an  arc  from  the  south  coast  of  Ko¬ 
rea  to  the  northwest  coast  of  Honshu.  Fang  and  Yang 
(1988),  Odamaki  (1989),  and  Teague  et  al.  (2001)  show 
that  tidal  currents  are  strong  in  many  portions  of  the  Ko¬ 
rea/Tsushima  Strait. 

Further  work  is  needed  to  obtain  sufficiently  accu¬ 
rate  tidal  height  and  current  predictions  for  the  Korea/ 
Tsushima  Strait.  Many  numerical  models  have  poor  reso¬ 
lution  in  the  Korea/Tsushima  Strait  and  therefore  cannot 
resolve  complex  tidal  structures  or  represent  the  interac¬ 
tion  of  the  tides  with  complex  topography.  Other  studies 
use  finer  resolution  models,  but  limit  the  number  of  tidal 
constituents  that  they  consider.  These  regional  modeling 
studies  often  force  their  models  at  the  open  ocean  bounda¬ 
ries  with  values  from  historical  co-tidal  charts  or  larger 
scale  models,  both  of  which  are  subject  to  error.  Obser¬ 
vational  studies  have  been  limited,  away  from  the  coast, 
to  short  time  series  of  tidal  current  data  in  which  some 
tidal  constituents  are  not  resolvable.  Moreover,  such  stud¬ 
ies  do  not  predict  tidal  currents  away  from  the  observa¬ 
tion  stations. 

Odamaki  (1989)  edited  the  tidal  charts  of  Ogura 
(1933)  based  on  a  large  observational  dataset  that  included 
data  from  many  short-duration  current  meters.  Odamaki 
(1989)  used  the  new  charts  produced  in  this  comprehen¬ 
sive  study  to  examine  the  tidal  dynamics  in  the  Korea/ 
Tsushima  Strait.  Recently,  Teague  et  al.  (2001)  have 


shown  that  the  Korea/Tsushima  tidal  charts  need  further 
adjustment.  They  analyzed  data  from  an  initial  deploy¬ 
ment  of  11  moorings,  each  equipped  with  an  acoustic 
Doppler  current  profiler  (ADCP)  and  a  pressure  gauge, 
along  two  lines  crossing  the  strait.  These  moorings  gath¬ 
ered  current  and  pressure  data  and  were  deployed  as  a 
part  of  the  U.S.  Naval  Research  Laboratory  (NRL)  “Link¬ 
ages  of  Asian  Marginal  Seas”  (LINKS)  program.  This 
program  will  be  described  in  more  detail  in  Section  2. 
Tidal  heights  were  obtained  by  harmonic  analysis  of  the 
pressure  time-series  and  then  compared  to  the  corrected 
chart  values  from  Odamaki  (1989).  Chart  amplitudes 
along  the  southern  line  of  moorings  are  about  10%  high 
for  the  semidiurnal  tides  and  10%  low  for  the  diurnal  tides. 
A  direct  comparison  between  the  tidal  current  coefficients 
obtained  from  a  LINKS  mooring  and  a  nearby  historical 
tidal  current  observation  show  a  40%  difference  in  M2 
amplitude  and  a  100%  difference  in  Oj  amplitude.  The 
most  likely  cause  for  these  differences  is  that  the  current 
observations  from  Odamaki  (1989)  were  of  too  short  du¬ 
ration  to  accurately  resolve  the  various  tidal  frequencies 
and  therefore  he  used  various  assumptions  to  account  for 
unresolved  constituents. 

This  paper  uses  the  tidal  characterizations  at  13  sta¬ 
tions  obtained  during  the  LINKS  study  and  other  data  with 
a  barotropic  data  assimilation  model  developed  by 
Griffin  and  Thompson  (1996).  The  model  produces  pre¬ 
dictions  of  tidal  currents  and  tidal  heights  by  using  the 
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Table  1.  Calculated  tidal  amplitudes  (amp.)  and  phases  (pha.)  in  cm  and  degrees  for  the  mooring  pressure  time  series  (Fig.  1). 
Nodal  corrections  were  applied  to  these  amplitudes  and  phases.  The  phases  are  the  Greenwich  phase  lags  according  to  the 
convention  given  by  Foreman  (1977).  The  estimated  amplitude  errors  (err.)  are  in  cm  and  were  calculated  as  described  in 
Subsection  3.2. 


NO  N2  N3  N4 


Oj  amp. 
Oj  err. 

O,  pha. 

P,  amp. 
P,  err. 

P,  pha. 
K,  amp. 
Kj  err. 
K,  pha. 
Lh  amP’ 

err. 

th  Pha- 
N2  amp. 
N2  err. 
N2  pha. 
M2  amp. 
M2  err. 
M2  pha. 
S2  amp. 
S2  err. 

S2  pha. 
K2  amp. 
K2  err. 

K2  pha. 


2.21 

0.10 

194.2 
0.58 
0.10 

259.5 

1.34 

0.10 

266.7 

1.02 

0.09 

334.4 

3.14 

0.09 

323.2 
16.28 

0.09 

333.0 

9.22 

0.09 

4.0 

2.69 

0.09 

353.0 


2.34 

0.06 

171.1 

0.26 

0.06 

197.6 

0.98 

0.06 

206.5 

0.92 

0.07 

348.8 

3.16 

0.07 

336.3 

17.08 

0.07 

344.0 

9.25 

0.07 

12.7 

2.48 

0.07 

3.2 


3.67 

0.07 

150.8 

0.81 

0.07 

159.6 
2.60 
0.07 

165.6 
1.05 
0.07 
1.8 
3.45 
0.07 

352.4 

18.28 

0.07 

357.8 

9.96 

0.08 

23.7 

2.75 

0.08 

14.6 


5.14 

0.08 

144.1 

1.35 

0.08 

156.7 
4.33 
0.08 

160.5 

1.16 

0.08 

10.7 
3.77 
0.08 
2.5 

19.86 

0.08 

7.2 

10.70 

0.08 

31.2 

2.93 

0.08 

22.0 


N5 

N6 

6.91 

9.12 

0.11 

0.10 

144.0 

144.9 

1.97 

2.66 

0.11 

0.10 

149.5 

162.9 

6.38 

8.43 

0.11 

0.10 

162.3 

166.7 

1.10 

1.39 

0.08 

0.09 

21.9 

24.0 

3.96 

4.59 

0.08 

0.09 

9.7 

19.9 

21.32 

23.43 

0.08 

0.09 

14.4 

23.9 

11.00 

12.29 

0.08 

0.09 

36.9 

45.0 

2.84 

3.42 

0.08 

0.09 

27.6 

35.8 

SI  S2 

14.28  13.34 

0.16  0.15 

34.7  44.1 

6.23  5.80 

0.17  0.15 

52.0  59.3 

19.88  18.41 

0.17  0.15 

58.5  65.8 

4.62  4.40 

0.16  0.14 

349.2  350.6 

17.84  16.53 

0.16  0.14 

352.5  352.1 

86.16  80.09 

0.16  0.14 

1.8  0.9 

40.14  37.52 

0.16  0.15 

26.1  25.2 

11.36  10.61 

0.16  0.15 

17.7  16.4 


S3  S4 

13.27  13.57 

0.11  0.12 

53.4  61.3 

5.64  5.65 

0.11  0.12 

67.4  74.7 

17.86  17.79 

0.11  0.12 

73.3  80.2 

4.19  4.08 

0.11  0.14 

351.6  353.0 

15.62  15.05 

0.11  0.14 

351.5  351.0 

76.09  73.87 

0.11  0.14 

0.1  359.0 

35.81  34.79 

0.11  0.14 

24.2  23.3 

10.15  9.79 

0.11  0.14 

15.1  14.2 


S5  S6 

14/78  16. 10 

0.13  0.14 

68.7  74.7 

5.89  6.32 

0.14  0.15 

79.8  88.7 

18.60  19.95 

0.14  0.15 

87.5  93.4 

3.40  4.04 

0.13  0.14 

357.0  356.5 

14.24  14.79 

0.13  0.15 

350.2  350.4 

72.11  72.85 

0.13  0.15 

357.4  358.0 

33.04  34.27 

0.13  0.15 

21.5  22.2 

8.73  9.71 

0.13  0.15 

12.9  12.6 


data  to  infer  optimal  open  ocean  boundary  conditions  for 
the  model.  This  paper  also  analyzes  the  accuracy  of  the 
predictions  and  estimates  the  accuracy  of  the  tide  model. 

2.  Tidal  Observations 

Figure  1  shows  the  locations  of  the  mooring  data  (tri¬ 
angles),  coastal  tide  station  data  (circles),  and  satellite 
altimeter  data  (diamonds)  used  in  this  work.  Thirteen  of 
the  sixteen  moorings  (N0-N6,  S1-S6,  Cl)  were  part  of 
the  LINKS  array.  Most  LINKS  sites  were  instrumented 
from  May  1999  through  March  2000  with  an  ADCP  and  a 
wave/tide  gauge,  housed  in  a  trawl-resistant  bottom 
mount.  The  ADCPs  measured  currents  throughout  the 
water  column  and  the  wave/tide  gauges  measured  pres¬ 
sure  at  the  bottom.  Due  to  mooring  failure,  instrument 
failure,  or  limited  numbers  of  instruments,  no  current  or 
pressure  data  are  available  from  NO  and  Cl  for  May 
through  October,  no  current  or  pressure  data  are  avail¬ 
able  from  S5  for  October  through  March,  and  no  pres¬ 
sure  data  are  available  from  Cl  and  N5  for  October 
through  March.  Further  details  of  this  mooring  effort  can 
be  found  in  Teague  et  ah  (2002).  Two  additional  bottom- 
mounted  ADCPs  were  deployed  during  the  time  of  the 


LINKS  moorings  by  the  Korea  Ocean  Research  and  De¬ 
velopment  Institute  (KORDI)  and  Pukyong  National  Uni¬ 
versity.  The  KORDI  mooring,  Kl,  was  deployed  in  the 
Jeju  (Cheju)  Strait  from  March  through  December  of 
1999.  The  Pukyong  National  University  mooring,  K2a  and 
K2b,  was  deployed  in  the  West  Channel  of  the  Korea/ 
Tsushima  Strait  from  June  through  December  of  1999. 
This  mooring  was  struck  by  fishing  gear  and  moved  at 
least  three  times.  The  data  from  the  KORDI  mooring  and 
the  data  from  the  first  two  positions  (June  through  Sep¬ 
tember,  October  through  November)  of  the  Pukyong  Na¬ 
tional  University  mooring  are  used  in  this  study.  The  data 
from  the  third  position  of  the  latter  mooring  are  not  used 
because  the  ADCP  failed  shortly  after  this  last  move. 

A  least-squares  harmonic  tidal  analysis  was  per¬ 
formed  on  all  current  and  pressure  data  measured  at  the 
mooring  sites.  Prior  to  this  analysis,  the  current  data  were 
vertically  averaged  to  obtain  estimates  of  the  barotropic 
currents.  The  tidal  analysis  was  done  using  the  Foreman 
(1977)  tidal  analysis  software  adapted  to  MATLAB  by 
Lentz,  Beardsley,  and  Pawlowicz.  The  analysis  by  Teague 
et  ah  (2001)  of  the  initial  deployment  (May  through  Oc¬ 
tober)  of  LINKS  moorings  found  that  the  eight  largest 
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tidal  constituents  in  the  pressure  data  were  O,,  P1(  K,,  fi2’ 
N2,  M2,  S2,  and  K2.  Excluding  low  frequency  and  shal¬ 
low  water  constituents,  these  are  also  the  eight  largest 
constituents  in  the  histogram  of  mean  tidal  amplitudes 
calculated  from  192  tide  gauges  in  the  East  China  Sea, 
the  Yellow  Sea,  and  the  Korea/Tsushima  Strait  by  Lefevre 
et  al.  (2000).  Therefore,  these  eight  constituents  were  used 
for  the  analysis  of  the  mooring  data  in  this  study.  The 
resulting  tide  parameters  are  presented  in  Tables  1-3. 
Records  from  K2a  and  K2b  were  too  short  to  resolve  the 


constituent  pairs  of  and  S2/K2,  and  therefore  P,  and 
K2  were  not  included  in  the  analysis  of  those  particular 
time  series. 

Both  the  ADCPs  and  wave/tide  gauges  in  this  study 
were  configured  so  that  each  measurement  was  an  aver¬ 
age  over  the  sampling  period  as  opposed  to  instantane¬ 
ous  samples.  The  Foreman  (1977)  tidal  analysis  package 
is  meant  to  analyze  instantaneous  samples  and  therefore 
corrections  are  required  to  prevent  bias  in  the  calculated 
amplitudes  and  phases.  Using  averages  instead  of  discrete 
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Table  3.  As  Table  2,  but  for  the  moorings  in  the  Jeju  Strait  and  vided  by  this  correction  factor,  which  is  very  nearly  equal 
in  the  West  Channel.  to  ^  for  these  data  because  At  values  were  either  15  or  30 
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values  causes  the  tidal  analysis  software  to  slightly  un¬ 
derestimate  tidal  amplitudes.  This  bias  is 


where  co  is  the  angular  frequency  of  the  tidal  constituent 
and  At  is  the  measurement  sampling  period.  The  ampli¬ 
tudes  calculated  by  the  tidal  analysis  software  were  di¬ 


minutes.  The  ADCPs  and  wave/tide  gauges  assign  the  start 
time  of  the  sampling  period  to  each  data  point.  To  obtain 
correct  phase  information  from  the  analysis,  each  data 
point  was  instead  assigned  the  time  of  the  middle  of  the 
sampling  period,  i.e.  adjusted  by  At/2. 

Differences  in  amplitudes  between  Tables  1  and  2 
and  the  tables  of  Teague  et  al.  (2001)  are  small  and  are 
mainly  because  this  study  uses  data  from  the  entire  dura¬ 
tion  of  the  mooring  deployments  as  opposed  to  the  dura¬ 
tion  of  the  initial  deployment.  However,  differences  in 
phases  are  not  small  because  the  time  adjustment  of  A//2 
required  by  the  interval  averaging  was  not  applied  in  the 
analysis  of  Teague  et  al.  (2001),  leading  to  errors  in  their 
tidal  phase  values. 

Further  tide  data  for  this  study  are  provided  by  a  tidal 
analysis  of  TOPEX/POSEIDON  (T/P)  data  and  from 
coastal  tide  station  data  (Fig.  1).  The  T/P  data  consist  of 
tidal  height  coefficients  calculated  from  selected  points 
on  tracks  107,  64,  18,  56,  and  102  using  the  method  de¬ 
scribed  in  Teague  et  al.  (2000).  This  calculation  estimates 
tidal  coefficients  using  aliased  frequencies  because  the 
T/P  altimeter  does  not  sample  quickly  enough  to  measure 
the  actual  tidal  frequencies.  Selection  was  done  to  reduce 
the  T/P  data  spatial  density  compared  to  that  of  the  moor¬ 
ing  data  and  varied  between  using  every  other  point  along 
the  track  to  using  every  fourth  point  along  the  track.  The 
coastal  tide  station  data  were  provided  as  constituent  co¬ 
efficients  and  are  from  the  International  Hydrographic 
Office  (IHO)  dataset.  Of  the  numerous  tidal  stations  lo¬ 
cated  in  the  Korea/Tsushima  Strait,  28  were  selected  for 
this  study  based  on  duration  of  their  measurements,  lack 
of  complexity  in  the  coastal  area  where  the  instrument 
was  placed,  and  overall  agreement  with  other  data  records. 

The  data  assimilation  model  (Section  3)  requires  time 
series  data.  Therefore,  the  Foreman  (1977)  tidal  analysis 
software  was  used  to  synthesize  time  series  from  the  co¬ 
efficients  obtained  by  the  analysis  of  the  mooring  data 
and  from  the  IHO  coefficients.  The  time  series  thus  gen¬ 
erated  are  for  a  period  of  41  days,  starting  in  May  of  1999, 
with  values  every  30  minutes.  Time  series  in  this  same 
form  were  synthesized  from  the  T/P  coefficients  using 
other  software  in  a  manner  consistent  with  their  coeffi¬ 
cient  format.  Some  of  the  IHO  stations  did  not  provide 
data  for  the  eight  constituents  considered  in  this  study, 
and  therefore  two  time  series  datasets  were  made.  The 
first  dataset  was  generated  at  all  data  sites,  but  using  only 
the  Oj,  Kj,  M2,  and  S2  coefficients.  The  second  dataset 
was  generated  from  all  eight  constituents  but  excluded 
those  sites  without  a  complete  set  of  constituent  data,  i.e. 
T3,  T7,  T8,  T13,  T14,  T16,  T17,  T20,  T23,  T24,  K2a,  and 
K2b  were  excluded.  These  two  datasets  provide  the  forc¬ 
ing  for  runs  of  the  data  assimilation  model. 
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3.  Barotropic  Data  Assimilation  Model 


3. 1  Model  dynamics 

A  numerical  model  was  used  to  assimilate  the  avail¬ 
able  data  and  provide  predictions  for  tidal  currents  and 
heights  throughout  the  Korea/Tsushima  Strait.  The  model 
uses  the  linearized  barotropic  shallow  water  equations; 


and 


du  dtl  ,  A 

^-  =  -gTL  +  ^'-T“’ 

dl  dx  h 


dv  dr]  ,  A 

-r-  =  STL-fr~Tv’ 

dt  dy  h 


dri  dhu  +  dhv 
dl  dy 


(2) 

(3) 


(4) 


where  the  terms  associated  with  the  incoming  wave  have 
been  combined  as  one  forcing  term,  F/h.  The  admittance 
coefficient  for  the  model  was  set  equal  to  -sjg  /  h ,  corre¬ 
sponding  to  the  admittance  of  a  plane  shallow  water  wave 
over  such  a  depth. 

The  open  ocean  boundaries  are  divided  into  a  series 
of  structure  functions.  Along  each  structure  function,  the 
term  F  from  Eq.  (7)  is  allowed  to  vary  in  time  but  is  forced 
to  remain  constant  along  its  segment  of  the  boundary.  A 
total  number  of  28  structure  functions  are  used,  spaced 
along  the  boundaries  such  that  the  bathymetry  is  relatively 
constant  for  each  structure  function.  They  varied  in  length 
from  5  to  75  km.  The  model  is  started  from  a  state  of  rest 
and  is  then  forced  by  the  28  time  series  of  F  along  the 
boundaries. 

3.2  Cost  function  and  weights 

A  cost  function  is  used  to  measure  how  well  a  par¬ 
ticular  model  run  successfully  predicts  the  observed  data 
set.  We  define  the  cost  function  to  be 


where  u  and  v  are  the  eastward  (x)  and  northward  (y) 
velocities,  respectively,  tj  is  the  surface  height,  t  is  time, 
g  is  the  gravitational  acceleration,  /  is  the  Coriolis  pa¬ 
rameter,  h  is  the  depth,  and  A  is  a  frictional  parameter. 

The  open  ocean  boundary  conditions  provide  the 
control  variables  for  the  model  as  now  described.  The 
boundary  condition  that  is  used  is  a  standard  Reid  and 
Bodine  (1968)  scheme.  This  scheme  can  be  interpreted 
in  the  following  way  for  the  model.  Consider  the  left  open 
ocean  boundary.  At  this  boundary,  the  normal  flow  can 
be  considered  to  be  a  superposition  of  an  incoming  (ex¬ 
ternally  forced)  gravity  wave  and  an  outgoing  (radiating) 
gravity  wave.  The  velocity  at  the  boundary,  u is  given 
by 

=  M;  +  m0,  (5) 

where  w(-  is  the  velocity  of  the  incoming  wave  and  u„  is 
the  velocity  of  the  outgoing  wave.  If  a  is  an  appropriate 
admittance  coefficient  for  these  gravity  waves,  then,  at 
the  left  boundary,  u0  equals  negative  a  multiplied  by  the 
surface  height  of  the  outgoing  wave,  r\0.  In  turn,  rja  is 
equal  to  the  total  surface  height  minus  the  surface  height 
due  to  the  incoming  wave,  77/-  Using  these  expressions, 
the  left  boundary  condition  becomes 

uL  =  Uj  -  a{r\L  -  Tjj)  (6) 

or 

uL=j-otr}L>  (?) 


- Pj,k)Wj(mhk  -pM),  (8) 

wherey  is  an  index  of  sites  and  types  of  data  to  be  used,  k 
is  the  temporal  index  of  the  data,  m  is  the  collection  of  all 
the  w,  v,  and  77  measurements,  p  is  the  corresponding  model 
predictions  for  these  measurements  produced  by  linear 
interpolations  from  the  model  space  and  time  grids,  and 
W  is  a  weighting  function.  The  best  fit  to  the  observed 
data  set  is  the  model  state  that  minimizes  the  cost  func¬ 
tion,  J,  assuming  that  model-data  differences  are  normally 
distributed  with  zero  means  and  W  is  determined  by  the 
inverse  of  the  error  variance  of  the  measurements 
(Thacker,  1987).  In  practice,  the  model  predictions, p,  are 
known  to  contain  bias  from  model  physics,  numerics,  and 
parameter  errors  (see  Subsection  4.4). 

The  cost  function  is  implicitly  a  function  of  the  F 
terms  from  Eq.  (7).  Choosing  values  of  F  at  each  struc¬ 
ture  function  for  each  control  time  completely  determines 
the  values  of  p  and  thus  J. .  Optimal  values  of  F  are  found 
by  minimizing  J ,  subject  to  the  strong  constraints  of  Eqs. 
(2)-(4)  bounded  by  Eq.  (7).  The  method  of  Lagrange 
multipliers  is  used  for  this  purpose,  producing  a  set  of 
adjoint  equations  that  the  multipliers  must  satisfy 
(Thacker,  1987).  The  system  of  equations  comprised  of 
the  discrete  forms  of  Eqs.  (2) — (4)  and  (7),  and  their 
adjoints  are  solved,  yielding  the  gradients  of  Eq.  (8)  with 
respect  to  each  structure  function  value  (F).  This  tech¬ 
nique  is  completely  described  in  Griffin  and  Thompson 
(1996)  and  summarized  in  Hallock  etal  (2003).  The  mini¬ 
mum  of  the  cost  function  is  then  found  using  this  gradi¬ 
ent  in  a  conjugate  gradient  method.  Values  of  F  are  started 
from  zero  and  then  adjusted  to  the  values  that  minimize 
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the  cost  function  through  this  process.  Convergence  was 
achieved  by  25  iterations  of  the  method. 

To  optimize  the  formulation  of  Eq.  (8),  the  weight¬ 
ing  of  the  measurements  must  be  optimized.  Assuming 
that  measurement  errors  are  uncorrelated  between  meas¬ 
urement  locations  and  are  stationary  in  time,  the  optimal 
weighting  for  Eq.  (8)  is 


where  crl  is  the  error  variance  of  the  measurements  and 
c  is  a  constant  term  for  all  values  of  j.  Unless  it  is  too 
small  or  too  large  for  the  computer’s  numerics  to  handle 
accurately,  the  value  of  c  does  not  affect  the  solution,  i.e. 
only  the  relative  weights  of  the  measurements  are  impor¬ 
tant. 

Estimates  of  <7^  ,  and  thus  WJt  were  made  using  the 
mooring  data  and  information  about  the  IHO  and  T/P  data. 
The  error  in  measuring  individual  velocity  or  pressure 
values  is  negligible  compared  to  the  error  in  separating 
the  tidal  oscillations  from  non-tidal  oscillations  in  the  data 
time  series.  Background  ocean  continuum  energy  causes 
uncertainty  in  the  least-squares  harmonic  analysis.  To 
estimate  the  level  of  this  energy,  spectra  were  calculated 
for  each  mooring  time  series,  high  energy  density  peaks 
were  set  to  the  median  energy  density  value  of  the  0.5- 
2.5  day""1  frequency  range,  and  a  straight  line  was  fitted 
to  the  adjusted  spectral  values  through  this  range.  Esti¬ 
mates  of  the  background  continuum  were  taken  as  values 
of  the  straight  line  at  Kj  and  M2  frequencies.  Errors  in 
the  least-squares  harmonic  analysis  for  the  diurnal  and 
semi-diurnal  constituents  were  calculated  in  the  standard 
way  using  a  variance  matrix  (calculated  from  the  time 
vectors  and  tidal  frequencies)  and  a  random  white-noise 
level  that  would  produce  the  same  energy  density  as  these 
calculated  Kj  or  M2  values.  Tables  1-3  present  the  esti¬ 
mated  amplitude  errors  for  the  mooring  data  using  this 
method. 

For  the  T/P  time  series,  the  variance  matrix  was  ob¬ 
tained  by  simulating  a  seven-year  time  vector  sampled  at 
the  T/P  repeat  frequency.  The  error  calculations  were  then 
done  using  this  variance  matrix  and  random  white-noise 
at  the  average  level  of  all  the  pressure  mooring  time-se¬ 
ries.  The  estimated  constituent  amplitude  errors  are  1.1— 
1.2  cm.  Potential  errors  due  to  non-tidal  energy  at  the 
aliased  frequencies  were  assumed  to  be  small  and  not 
considered.  The  periods  of  these  aliased  frequencies  in 
days  are  45.71  for  0„  88.89  for  P„  173.19  for  K1?  20.31 
for  p2, 49.53  for  N2, 62. 1 1  for  M2, 58.74  for  S2,  and  86.60 
for  K2. 

For  the  IHO  coefficients,  time  vectors  were  simu¬ 
lated  with  1  hour  sampling  for  the  duration  given  in  the 


IHO  database.  The  error  calculations  were  done  consid¬ 
ering  all  the  diurnal  and  semi-diurnal  constituents  present 
in  the  IHO  record  to  form  the  variance  matrix  and  using 
random  white-noise  at  the  average  level  of  all  the  pres¬ 
sure  mooring  time-series.  At  some  short  duration  stations, 
the  main  eight  constituents  plus  additional  constituents 
were  included  in  the  analysis  despite  violation  of  the 
Rayleigh  criteria  by  constituents  close  in  frequency.  The 
standard  error  analysis  using  the  variance  matrix  accounts 
for  the  magnification  of  background  error  caused  by  this 
inclusion.  Alternatively,  at  other  short  duration  stations, 
some  of  the  main  eight  constituents  were  not  included  in 
the  analysis.  At  these  sites,  the  non-analyzed  or  missing 
constituent  of  the  Rayleigh  pair  will  interfere  with  the 
analyzed  constituent,  leading  to  errors  proportional  to  the 
amplitude  of  the  missing  constituent.  For  the  error  calcu¬ 
lations,  the  unknown  amplitude  of  the  missing  constitu¬ 
ent  was  taken  as  the  model  prediction  from  the  eight-con¬ 
stituent  run.  This  same  procedure  was  done  to  account 
for  errors  at  moorings  K2a  and  K2b  caused  by  the  miss¬ 
ing  Pj  and  K2  constituents. 

Errors  caused  by  limited  precision  in  the  reported 
IHO  coefficients  were  also  estimated  using  the  model 
predictions  from  the  eight-constituent  run  at  the  data  sites. 
At  several  sites  the  precision  error  was  larger  than  the 
analysis  error  and  therefore  was  used  in  place  of  the  analy¬ 
sis  error.  The  final  estimated  error  for  the  IHO  coeffi¬ 
cients  varied  widely  from  site  to  site  with  a  maximum  at 
station  T7  of  6.00  cm  for  S2  and  a  minimum  at  station 
T28  of  0.06  cm  for  each  constituent.  The  former  error  is 
large  because  of  the  lack  of  inclusion  of  K2  in  the  analy¬ 
sis  of  this  3  month  record  and  the  latter  error  is  small 
because  the  station  at  Hamada  was  of  10  year  duration 
and  had  high  precision  in  the  reported  results. 

Values  for  G2m  were  calculated  from  the  estimated 
amplitude  errors  of  the  constituents  and  then  used  in  Eq. 
(9)  to  calculate  separate  Wj  for  the  four  and  eight  con¬ 
stituent  runs.  Numerous  preliminary  runs  were  done, 
changing  various  weights  to  test  the  model’s  sensitivity 
to  them.  For  most  of  the  weight  adjustments,  changes  to 
the  final  solution  were  small.  However,  in  some  regions 
of  the  model,  such  as  south  of  Jeju  Island,  changing  the 
weights  can  significantly  alter  model  predictions. 

3.3  Model  bathymetry  and  friction  parameter 

The  model  domain  is  an  Arakawa  C  grid  with  a  2.5 
km  grid  spacing.  The  bathymetry  field,  h,  is  defined  every 
2.5  km,  while  «,  v,  and  7]  are  defined  every  5  km.  Several 
different  bathymetries  were  used  in  preliminary  runs.  The 
model  was  first  run  at  a  lower  resolution  (5  km)  using 
both  the  publicly  available  U.S.  Navy  DBDB-V 
bathymetry  database  (Naval  Oceanographic  Office,  1997) 
and  a  bathymetry  field  (SKKU  dataset)  produced  from  a 
1 -minute  resolution  data  set  available  from  the  Labora- 
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tory  for  Coastal  and  Ocean  Dynamics  Studies,  Sung  Kyun 
Kwan  University  (Choi,  1999).  The  value  of  the  final  cost 
function  was  52%  higher  for  the  DBDB-V  dataset,  and 
therefore  of  these  two  datasets,  only  the  SKKU  dataset 
was  used  at  the  2.5  km  resolution  runs.  The  high  resolu¬ 
tion  runs  with  the  SKKU  dataset  produced  good  results 
but  some  position-dependent  differences  remained  be¬ 
tween  the  model  predictions  and  the  data,  seemingly  re¬ 
lated  to  bathymetry.  The  SKKU  dataset  does  not  have  deep 
enough  values  for  the  deep  trough  that  exists  in  the  West 
Channel  between  Tsushima  and  Korea.  For  example,  ob¬ 
served  depths  for  the  two  positions  of  K2  were  228  and 
156  m  respectively,  but  the  SKKU  dataset  gives  depths 
of  only  161  and  122  m.  The  SKKU  dataset  also  appears 
to  be  inaccurate  near  mooring  site  SI.  Further,  the  SKKU 
dataset  does  not  include  several  of  the  larger  islands  along 
the  south  coast  of  Korea. 

Although  use  of  the  SKKU  dataset  produces  good 
results  that  are  significantly  better  than  using  the  DBDB- 
V  dataset,  this  study  used  a  third  bathymetric  database 
for  the  final  model  runs.  This  final  bathymetry  database 
was  formed  by  gridding  (Smith  and  Wessel,  1990) 
bathymetry  soundings  obtained  from  the  National  Imagery 
and  Mapping  Agency  (NIMA)  Digital  Nautical  Charts 
(DNC®).  Land  masking  was  done  using  the  World  Vector 
Shoreline  land  mapping  from  the  Generic  Mapping  Tool 
(Wessel  and  Smith,  1998)  at  high  resolution  and  fitting 
the  results  to  the  2.5  km  grid  of  this  study.  The  resulting 
final  bathymetric  database  contains  finer  resolution  fea¬ 
tures  than  either  the  DBDB-V  or  the  SKKU  datasets,  bet¬ 
ter  represents  the  true  depths  in  the  West  Channel  trough, 
and  produces  lower  minimized  cost  function  values  than 
runs  made  with  the  SKKU  dataset. 

The  value  of  A  in  Eqs.  (2)  and  (3)  was  varied  be¬ 
tween  2.3  x  10"4  and  2.3  x  10“2  m/s  for  different  model 
runs.  The  value  of  the  cost  function  decreases  and  then 
increases  as  the  friction  parameter  is  lowered,  and  a  value 
of  A  =  6.8  x  10"4  m/s  was  chosen  for  the  final  model  runs 
because  this  parameter  value  produced  the  lowest  cost 
function  value.  The  functional  relationship  between  cost 
function  value  and  frictional  parameter  is  nearly  flat  near 
this  point  (for  example,  A  =  4.6  x  10-4  m/s  produces  nearly 
the  same  cost  function  value)  but  rises  steeply  for  A  val¬ 
ues  greater  than  1.1  x  10“3  m/s.  The  magnitude  of  the 
slope  also  increases  for  A  values  lower  than  4.6  x  10-4 
m/s  but  this  trend  could  not  be  fully  evaluated  because 
the  numerics  of  the  model  failed  for  friction  parameters 
lower  than  2.3  x  10-4  m/s.  Kang  et  al  (1995)  estimated 
linearized  bottom  friction  coefficients  in  different  sea¬ 
sons  for  the  Korea/Tsushima  Strait.  Their  early  spring 
value  (converted  to  the  form  used  in  this  paper)  of  1.3  x 
10~3  m/s  and  their  summer  value  of  1.8  x  10-3  m/s  are 
both  higher  than  our  optimized  value. 

In  summary,  this  data  assimilation  technique  attempts 


Fig.  2.  Rms  differences  between  the  observation  time  series 
and  the  model  prediction  time  series.  Open  circles  indicate 
rms  differences  from  the  four-tidal  constituent  run.  Stars 
indicate  rms  differences  from  the  eight-constituent  model 
run.  Stars  are  not  plotted  for  sites  that  lack  data  for  all  eight 
constituents.  See  Fig.  1  for  locations  of  the  data  sites. 


to  match  observed  values  of  «,  v,  and  rj  by  adjusting  the 
forcing  along  the  model  boundary.  The  model  is  con¬ 
strained  to  follow  the  discrete  representation  of  the  dy¬ 
namical  equations  at  all  the  internal  grid  points.  Further 
details  of  this  model  and  data  assimilation  technique  are 
presented  in  Griffin  and  Thompson  (1996). 

4.  Results  and  Error  Analysis 

4. 1  Model/data  comparisons 

How  well  does  the  data  assimilation  model  predict 
the  observations?  This  question  could  be  answered  by 
examining  the  minimized  cost  function  (Eq.  (8)).  How¬ 
ever,  it  is  more  revealing  to  examine  the  results  of  a  root 
mean  square  (rms)  difference  calculation,  as  shown  in  Fig. 
2.  This  figure  shows  the  rms  of  the  observation  time  se¬ 
ries  minus  the  model  predicted  time  series  at  all  data  sites 
for  both  the  four-constituent  run  (circles)  and  the  eight- 
constituent  run  (stars).  The  rms  differences  for  both  runs 
are  very  similar. 

Rms  differences  between  the  measurements  and  pre¬ 
dictions  of  either  component  of  velocity  are  below  4.3 
cm/s  at  all  sites  and  below  2.4  cm/s  at  all  sites  except 
K2a  and  K2b.  Although  located  close  to  each  other,  the 
tidal  parameters  measured  at  Cl,  K2a,  and  K2b  all  had 
different  values  (Table  3).  The  differences  between  the 
measurements  at  Cl  and  K2a  are  exceptionally  large. 
Measurements  made  at  K2a  and  K2b  have  much  higher 
Kj  and  S2  errors  than  measurements  at  other  sites,  but  the 
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Fig.  3.  Model  predicted  co-tidal  (solid)  and  co-range  (dashed) 
lines  for  the  M2  tide.  The  co-range  lines  are  in  meters  with 
increments  of  0.1  m.  The  co-tidal  lines  are  in  degrees  with 
increments  of  15  degrees.  The  Greenwich  phase  conven¬ 
tion  used  in  this  figure  is  described  by  Foreman  (1977). 


Fig.  5.  As  Fig.  3,  but  for  the  Kj  tide.  The  co-range  lines  are  in 
increments  of  0.025  m.  The  co-tidal  lines  are  in  increments 
of  45  degrees.  The  gray  circle  indicates  the  position  of  an 
amphidromic  point 


Fig.  4.  As  Fig.  3,  but  for  the  S2  tide.  The  co-range  lines  are  in 
increments  of  0.05  m.  The  co-tidal  lines  are  in  increments 
of  15  degrees. 


Fig.  6.  As  Fig.  3,  but  for  the  Ox  tide.  The  co-range  lines  are  in 
increments  of  0.025  m.  The  co-tidal  lines  are  in  increments 
of  45  degrees.  The  gray  circle  indicates  the  position  of  an 
amphidromic  point. 


differences  between  the  measured  values  at  Cl  and  K2a 
exceed  these  errors.  These  measured  small-scale  differ¬ 
ences  may  be  caused  by  bathymetric  effects  at  the  deep 
trough  where  these  moorings  were  deployed.  Alterna¬ 
tively,  it  is  also  possible  that  an  unknown  source  of  error 
has  contaminated  the  measurements  at  K2a.  Regardless, 
because  the  estimated  errors  for  Cl  are  lower  than  those 
for  K2a  and  K2b,  the  data  assimilation  model  gave  pref¬ 
erence  to  the  Cl  data. 

The  average  value  of  the  velocity  component  rms 
differences  is  1.2  cm/s  for  the  eight  constituent  run  and 

1.4  cm/s  for  the  four  constituent  run.  The  average  values 
of  the  rms  differences  for  surface  height  are  0.9  cm  for 
the  LINKS  mooring  sites,  3.2  cm  for  the  coastal  tide  sta¬ 
tions,  and  4.3  cm  for  the  T/P  sites  using  the  predictions 
from  the  eight-constituent  run.  For  the  four-constituent 
run,  the  values  are  0.9  cm  for  the  LINKS  mooring  sites, 

3.4  cm  for  the  coastal  tide  stations,  and  3.7  cm  for  the  T/ 
P  sites.  Highest  differences  occur  along  the  west  coast  of 
the  Gotou  Islands,  off  the  southeast  coast  of  Jeju  Island, 


in  the  West  Channel,  and  along  the  edge  of  the  southern 
model  boundary. 

Differences  between  observation  and  modeled  time 
series  are  small  compared  to  the  total  tidal  fluctuations 
in  the  Korea/Tsushima  Strait.  The  average  standard  de¬ 
viations  for  the  observed  eight-constituent  tidal  time  se¬ 
ries  are  13.1  cm/s  for  the  velocity  time  series,  38.9  cm 
for  the  surface  height  time  series  at  the  LINKS  moor¬ 
ings,  49. 1  cm  for  the  surface  height  time  series  at  the  tidal 
stations,  and  42.5  cm  for  the  surface  height  time  series  at 
the  T/P  sites.  The  normalized  rms  difference  for  the  ve¬ 
locity  is  less  than  0.15  for  all  time  series  except  for  the  v 
time  series  at  SI,  S2,  Kl,  and  K2a  and  the  four-constitu¬ 
ent  u  time  series  at  S6  and  K2b.  Variances  of  the  observed 
v  time  series  at  SI,  S2,  and  Kl  are  small.  The  normalized 
rms  difference  for  surface  height  is  less  than  0.15  for  all 
time  series  except  for  time  series  at  sites  A2-A3  and  at 
sites  A23-A29.  Variances  of  the  observed  time  series  at 
A23-A29  are  small. 
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Fig.  7.  Model  predicted  tidal  current  ellipses  for  the  M2  tide. 
Linear  interpolation  was  used  over  the  model’s  C  grid  in 
order  to  collocate  u  and  v  predictions.  For  display  purposes, 
every  fourth  ellipse  is  drawn.  Three  ellipses  in  the  north¬ 
west  Jeju  Strait  and  one  ellipse  in  the  narrows  between  the 
Gotou  Islands  and  Kyushu  are  drawn  at  half  scale  with  a 
dashed  line.  Dots  are  drawn  to  indicate  the  tip  of  the  tidal 
current  vector  at  the  time  when  the  astronomical  argument 
is  zero  (Foreman,  1978).  Stars  are  drawn  in  place  of  dots 
for  the  half  scale  ellipses. 


Fig.  8.  As  Fig.  7,  but  for  the  Kt  tide.  All  ellipses  are  drawn 
with  the  same  scale. 


4.2  Model  predictions 

Figures  3-12  present  the  model  tidal  predictions 
throughout  the  model  domain.  The  fields  were  made  by 
harmonically  analyzing  the  model’s  time  domain  output. 
Nodal  corrections  have  been  applied  to  all  amplitudes  and 
phases.  Fields  for  M2,  S2,  K,„  and  O,  were  generated  from 
the  four-constituent  model  run;  those  for  N2,  Pj,  K2,  and 
ji2  were  generated  from  the  eight-constituent  model  run. 
The  Rayleigh  criterion  is  violated  for  the  Kj/Pi  and  the 
S2/K2  pairs  in  the  harmonic  analysis  of  this  latter  run  be¬ 
cause  of  the  41-day  duration  of  the  model  time  series. 
However  this  does  not  increase  the  error  for  these  con¬ 
stituents  because  the  background,  i.e.  non-tidal,  energy 
from  the  model  is  negligible.  Tidal  amphidromes  are 
marked  in  these  figures  if  the  amphidrome  is  completely 
enclosed  by  the  model  domain. 


Fig.  9.  As  Fig.  7,  but  for  the  S2  tide.  All  ellipses  are  drawn  with 
the  same  scale. 


Fig.  10.  As  Fig.  7,  but  for  the  0{  tide.  All  ellipses  are  drawn 
with  the  same  scale. 


4.3  Model  prediction  accuracy 

The  model  provides  tidal  predictions  at  all  grid  points 
as  well  as  at  the  observation  locations.  How  well  does 
the  model  represent  the  true  barotropic  tide  throughout 
its  domain?  For  the  purposes  of  further  error  analysis  it 
is  convenient  to  express  the  data  assimilation  model  in  a 
matrix  form,  as  was  done  by  Dowd  and  Thompson  (1996) 
and  Thompson  et  al.  (2000).  Thompson  et  ah  (2000) 
showed  that  for  one  tidal  constituent,  a  form  of  the  model 
can  be  written  as 


^  =  (10) 

where  pD  is  a  column  vector  containing  sine  and  cosine 
amplitude  coefficients  from  the  harmonically  analyzed 
model  output  at  all  observation  locations,  RD  is  a  response 
matrix,  and  ft  is  a  column  vector  containing  sine  and  co¬ 
sine  amplitude  coefficients  for  all  the  structure  functions 
from  the  harmonically  analyzed  forcing  time  series  (F 
from  Eq.  (7)).  Following  the  same  methods  presented  in 
the  appendices  of  Dowd  and  Thompson  (1996),  another 
form  of  the  model  can  be  written  as 


A =  (11) 
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Fig.  11.  As  Fig.  3,  but  for  the  N2,  Pj,  K2,  and  fa  tides.  The  co-range  lines  are  in  increments  of  0.025  m,  0.020  m,  0.020  m,  and 
0.015  m  respectively.  The  co-tidal  lines  are  in  increments  of  30  degrees  for  N2,  K2,  and  fa  and  in  increments  of  45  degrees  for 
Pj.  The  gray  circle  indicates  the  position  of  the  Pj  amphidromic  point. 
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Fig.  12.  As  Fig.  7,  but  for  the  N2,  P,,  K2,  and  fa  tides. 
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ellipses  are  all  at  5  times  the  scale  used  in  Figs.  8-10. 


where  jJ.G  is  a  column  vector  containing  sine  and  cosine 
amplitude  coefficients  from  the  harmonically  analyzed 
model  output  at  all  grid  points,  RG  is  a  different  response 
matrix,  and  /3  is  as  defined  above.  Equations  (10)  and  (11) 
represent  the  forward  model  for  a  particular  tidal  con¬ 
stituent.  The  actual  model  is  run  in  the  time  domain  as 
described  in  Section  3.  Solutions  for  RD  and  RG  were  ob¬ 
tained  from  time  domain  runs  of  the  forward  model  by 
following  the  alternative  method  described  at  the  end  of 
appendix  A  in  Dowd  and  Thompson  (1996).  This  method 


forces  the  model  with  a  unit  amplitude  sine  or  cosine  time 
series  at  a  single  structure  function  and  zero  at  all  other 
structure  functions.  In  matrix  form,  one  value  of  (5  has 
been  set  to  one  and  all  other  values  set  to  zero.  The  har¬ 
monically  analyzed  output  from  such  a  run  yields  a  par¬ 
ticular  column  of  RD  and  a  particular  column  of  RG .  Fur¬ 
ther  columns  of  RD  and  RG  were  obtained  from  additional 
forward  model  runs  by  changing  the  structure  function 
used  for  forcing. 

Thompson  et  al.  (2000)  state  the  matrix  form  of  their 
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unweighted  cost  function  and  then  calculate  the  optimal 
choice  for  /?  and  the  expected  variance  of  this  optimal 
choice  using  standard  linear  regression  techniques.  Us¬ 
ing  our  weighted  cost  function  (Eq.  (8)),  the  matrix  equa¬ 
tions  for  the  cost  function,  the  optimal  choice  for  /J,  and 
the  expected  variance  of  this  optimal  choice  are 

J  =  (m-tiD)'w(m-nD),  (12) 

p0=(R'DwRDYlR'Dwm,  (13) 

and 

var(j80)  =  (R'DwRD)~l  R'Dwvar(m)wRD(R'DwRDf\  (14) 

where  m  is  a  column  vector  containing  the  sine  and  co¬ 

sine  amplitude  coefficients  for  the  harmonically  analyzed 
observation  time  series,  fi0  is  the  optimal  choice  for  /?, 
and  w  is  a  diagonal  matrix  containing  weighting  values 
derived  from  the  values  of  W  in  Eq.  (8).  Primes  indicate 
the  transpose  operation.  The  notation  var(-)  will  be  used 
throughout  the  paper  to  indicate  the  ensemble  variance 
of  a  quantity. 

For  the  model  used  in  this  study,  there  are  particular 
values  of  /J  such  that  fiD  in  Eq.  (10)  is  approximately  zero. 
The  amplitudes  and  the  error  variances  of  the  amplitudes 
for  such  singular  vectors  will  be  ill  determined  in  Eqs. 
(13)  and  (14)  and  reach  high  values.  Nevertheless,  final 
amplitudes  of  these  ill-determined  vectors  in  the  time 
domain  model  remain  low  because  of  the  initialization  of 
the  vectors  to  zero  and  the  conjugate  gradient  methodol¬ 
ogy.  Therefore,  the  matrix  model  analog  requires  a  spe¬ 
cial  treatment  of  these  vectors  to  accurately  reflect  the 
time  domain  model  and  to  prevent  unrealistic  forcings 
from  affecting  the  error  estimates. 

To  accomplish  this,  first  the  singular  vectors  of  RD 
were  found  by  a  singular  value  decomposition, 

Rd=UKV.  (15) 

V  is  a  unitary  matrix  with  orthonormal  columns  that  are 
the  singular  vectors  of  RD.  These  vectors  are  ranked  by 
row  according  to  the  amount  of  response  produced  at  the 
data  sites  by  using  them  as  /?  in  Eq.  (10).  The  response 
levels  are  quantified  by  the  diagonal  of  A  and  the  struc¬ 
ture  of  the  responses  are  given  by  the  columns  of  U.  Con¬ 
verting  the  problem  to  the  solution  space  spanned  by  the 
singular  vectors  gives 

R*=RdV  =  UA.  (16) 

Now  the  solution  space  can  be  divided  into  the  determined 
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Fig.  13.  Estimated  error  in  the  predictions  of  M2  rj  tidal  ampli¬ 
tudes  from  the  deterministic  part  of  the  solution  space  (see 
Subsection  4.3).  Contours  are  drawn  at  0.005  m  increments. 
The  0.005  and  0.01  m  contours  are  drawn  as  bold  lines 


solution  space  and  the  null  space.  For  this  study,  the  di¬ 
viding  criteria  classified  a  singular  vector  in  the  null  space 
if  all  responses  to  the  vector  at  the  data  sites  were  less 
than  1%  of  the  vector  amplitude,  i.e.  the  corresponding 
column  of  RD*  had  values  all  less  than  0.01.  For  M2,  18 
of  the  singular  vectors  were  classified  in  the  null  space 
and  38  vectors  were  classified  in  the  determined  space. 
After  this  separation,  the  null  space  can  be  truncated  from 
the  solution  by  defining  a  matrix  RD*  that  is  made  up  of 
the  first  38  columns  of  and  a  matrix  RD f  that  is  made 

up  of  the  last  18  columns  of  RD *.  In  this  new  space,  the 
variances  of  the  determined  singular  vectors,  /}*,  are  found 
by  applying  Eq.  (14)  using  RD *  instead  of  RD.  The 
variances  of  the  null  singular  vectors,  are  set  to  one 
and  their  co-variances  are  set  to  zero.  Large  values  for 
the  amplitudes  of  these  vectors  would  either  produce  de¬ 
tectable  changes  at  the  data  sites  or  lead  to  unrealisti¬ 
cally  large  tidal  fluctuations  elsewhere  in  the  model  do¬ 
main  and  therefore  cannot  occur  in  nature. 

In  Eq.  (11),  the  tidal  predictions  at  the  grid  points 
are  determined  by  a  constant  term  multiplied  by  j80.  There¬ 
fore,  in  the  new  space,  the  expected  variance  matrices 
for  these  predictions  are 


var 


=  wRp  j  R-d  wveLT(m)wR*D^R*D  V*  R'a 


(17) 


and 


var(nl)  =  RGvUv(p')v''R'G,  (18) 
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Fig.  14.  Estimated  error  in  the  predictions  of  M2  77  tidal  ampli¬ 
tudes  caused  by  the  undetermined  singular  vectors  (see 
Subsection  4.3).  The  0.005  and  0.01  m  contours  are  drawn 
as  bold  lines.  The  0.02,  0.03,  0.04,  and  0.05  m  contours  are 
drawn  with  normal  thickness.  Additional  contours  are  drawn 
at  0.1  m  increments. 


where  F*  and  F1*  are  the  first  38  columns  of  F and  the  last 
18  columns  of  F  respectively.  Figures  13  and  14  present 
the  predicted  M2  rj  errors  for  the  four-constituent  run  re¬ 
sulting  from  the  use  of  Eqs.  (17)  and  (18).  The  predicted 
errors  for  most  of  the  model  domain  are  less  than  0.5  cm. 
The  model  null  space  (Fig.  14)  occurs  along  the  far  south¬ 
western  coast  of  Korea,  south  of  Jeju  Island,  between  the 
Gotou  Islands  and  Kyushu,  and  along  the  coast  of  Honshu 
north  of  35.2°.  The  errors  shown  in  this  figure  are  true 
only  if  the  amplitudes  of  the  null  space  vectors  can  in 
actuality  reach  values  of  one. 

4.4  Model  accuracy 

The  theory  in  Subsection  4.3  has  assumed  that  the 
model  physics,  numerics,  and  parameters  are  true.  Under 
these  assumptions  the  error  in  the  model  over  most  of  the 
domain  is  very  small.  How  large  are  the  errors  in  the 
model  due  to  imperfect  model  physics,  numerics,  and 
parameters?  Figure  2  suggests  that  these  errors  may  be 
larger  than  the  errors  calculated  from  Eq.  (17). 

For  ease  of  notation,  let 

S  =  Rd”-  (19) 

Using  this  expression  and  combining  Eqs.  (10)  and  (13), 
the  measured  M2  amplitude  values  minus  their  model 
counterparts  are  then 

m  —  liD~  m  -  Sm  =  (7  -  S)m.  (20) 

In  this  equation,  I  is  the  identity  matrix.  Using  the  trun¬ 
cated  form  of  Rd  in  Eq.  (19)  allows  the  accurate  evalua¬ 
tion  of  the  matrix  inverse  and  does  not  drastically  affect 


Fig.  15.  Differences  between  the  observed  and  modeled  M2 
cosine  (circles)  and  sine  (pluses)  amplitudes  at  all  the  data 
sites.  Dashed  lines  are  drawn  ±  one  standard  deviation  from 
the  expected  value  (zero)  of  these  differences.  The  posi¬ 
tions  of  these  lines  were  determined  from  Eq.  (22). 


the  results  because  the  null  space  vectors  have  small  pro¬ 
jections  to  the  data  sites.  If  m  is  composed  of  random 
noise  plus  a  signal  that  satisfies  the  model  physics, 
numerics,  and  parameters,  then  the  expected  value,  E(-), 
and  variance  of  Eq.  (20)  are 

E(m-fiD)  =  0  (21) 

and 

var [m  -  flD)  -  var (m)  -Svar(m)  -  var (tn)S'  +  Svzr(m)S' . 

(22) 

In  Eq.  (22),  the  second  and  third  terms  are  equal  because 
S  is  symmetric  and  var(w)  is  diagonal.  The  magnitude  of 
the  fourth  term  is  also  equal  to  the  magnitude  of  these 
terms  if  the  weighting  in  the  model  is  optimal,  var (m) 
represents  the  random  error  in  the  data  at  each  site  and 
-Svar (m)  represents  over  fitting  of  the  model.  The  com¬ 
bined  term  of  -var (m)^  +  represents  the  affect 

of  weighting  some  sites  too  strongly  and  others  too 
weakly. 

Figure  15  shows  the  difference  between  the  measured 
M2  sine  and  cosine  amplitudes  and  the  predicted  M2  am¬ 
plitudes  for  the  four-constituent  run.  Values  ±  one  stand¬ 
ard  deviation  from  the  expected  value  (zero)  of  these  dif¬ 
ferences,  calculated  from  Eq.  (22),  are  plotted  as  dashed 
lines.  If  the  theory  presented  by  Eqs.  (21)  and  (22)  is  cor¬ 
rect  then  the  chance  is  small  that  the  differences  in  Fig. 
15  would  be  far  above  or  below  one  standard  deviation. 
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Fig.  16.  Estimated  bias  in  the  predictions  of  M2  t]  tidal  ampli¬ 
tudes.  Contours  are  drawn  at  0.01  m  increments.  Negative 
values  (dashed  contours)  indicate  that  the  predicted  model 
amplitudes  are  too  low.  These  error  values  were  derived 
from  Eq.  (25). 


Fig.  17.  The  total  estimated  rms  error  for  the  M2  r/  tidal  ampli¬ 
tudes.  These  errors  are  equal  to  the  square  root  of  the  sum 
of  the  squares  of  the  errors  in  Figs.  13,  14,  and  16.  The 
0.005  m  contour  is  drawn  as  a  dashed  line  and  the  0.01  m 
contour  is  drawn  as  a  bold  line.  The  0.02,  0.03,  0.04,  and 
0.05  m  contours  are  drawn  with  normal  thickness.  Addi¬ 
tional  contours  are  drawn  at  0.1  m  increments. 


Statistically,  this  theory  can  be  rejected  at  the  99%  confi¬ 
dence  level  at  all  the  velocity  sites  and  at  54%  of  the  rj 
sites.  Most  of  the  sites  where  the  theory  cannot  be  re¬ 
jected  based  on  the  evidence  are  T/P  measurement  loca¬ 
tions. 

There  are  two  possible  explanations  for  the  failure 
of  the  theory.  The  first  is  that  the  error  variances  esti¬ 
mates,  var(w),  are  inaccurate  and  the  standard  deviation 
lines  in  Fig.  15  should  be  farther  away  from  zero.  The 
second  is  that  m  is  composed  of  more  than  random  noise 
plus  tidal  signals  that  follow  the  linear  model  physics  and 
model  parameters,  i.e.  the  model  physics,  numerics,  and/ 
or  parameters  contain  error.  The  second  explanation  is 
more  likely  because  of  the  large  degree  of  departure  from 
the  error  theory  at  many  sites  where  the  error  variance 
estimates  should  be  very  accurate,  and  because  of  the 
many  possible  sources  of  error  including  bathymetry  er¬ 
rors,  error  in  the  friction  formulation,  error  in  the  friction 
coefficient,  resolution  of  the  structure  functions,  non-lin¬ 
ear  tide  interactions,  etc.  Therefore,  this  study  will  as¬ 
sume  that  the  error  variance  estimates  are  accurate  and 
that  the  larger  errors  are  due  to  errors  in  the  model,  or 
errors  in  the  model  parameters. 

To  account  for  this,  let  m  now  be  composed  of  a  sig¬ 
nal  that  satisfies  the  model,  plus  a  true  tide  signal,  NDl 
that  is  not  represented  in  the  model,  plus  random  noise. 
With  this  new  theory,  the  variance  of  ( m  -  / iD )  is  still 
given  by  Eq.  (22)  but  Eq.  (21)  becomes 

E(m  -  }iD)  =  (/-  S)Nd.  (23) 

Although  the  data  provide  only  a  single  realization  of 
(m  -  jJ,D)t  this  single  realization  can  be  used  as  an  accu¬ 
rate  estimate  of  the  expected  value  because  the  variance 
of  this  quantity  is  low  (Fig.  15).  The  matrix  (/  -  S)  is 
always  singular  because  of  the  structure  of  S  and  there¬ 


fore  Eq.  (23)  cannot  be  simply  inverted  to  solve  for  the 
unknown  values  of  ND.  However,  an  approximation  to  ND 
was  found  by  solving  the  related  equation 

(/- S)\m  -  fiD)  =  (/-£)'( I-S)ND  (24) 

using  the  conjugate  gradient  method  (Golub  and  Van 
Loan,  1996). 

The  bias  of  jlG  is 

bias(/zc)  =  tfcr  (<><>)  Rp'wND  -  Ng,  (25) 

where  NG  is  the  true  tide  signal  at  the  model  grid  points 
that  is  not  represented  by  the  model.  NG  was  estimated 
by  gridding  (Smith  and  Wessel,  1990)  the  values  of  ND 
and  then  Eq.  (25)  was  used  to  estimate  the  bias  in  jaG. 
Figure  16  shows  the  M2  amplitude  bias  estimate  derived 
from  the  results  of  Eq.  (25).  Negative  values  (dashed  con¬ 
tours)  indicate  that  the  predicted  model  amplitudes  are 
too  low.  In  general  the  M2  bias  field  is  between  -2  and  1 
cm,  but  it  reaches  values  less  than  -3  cm  to  the  west  of 
Busan,  around  Tsushima  Island,  to  the  east  of  Jeju  Island, 
and  along  the  coast  of  Kyushu.  Values  higher  than  2  cm 
occur  directly  to  the  west  of  the  Gotou  Islands  and  in  the 
null  space  south  and  north  of  Jeju  Island. 

Figure  17  shows  the  estimated  rms  error  for  the  M2 
amplitude  field.  This  map  combines  the  errors  from  Figs. 
13,  14,  and  16,  but  the  errors  from  Fig.  16  dominate  in 
most  locations.  For  broad  regions  near  the  north  and  south 
LINKS  mooring  lines  and  in  the  Jeju  Strait  the  rms  errors 
are  less  than  0.5  cm.  Elsewhere,  outside  the  null  space, 
errors  range  from  1  to  4  cm. 
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5.  Discussion 

Figures  3-12  show  the  complex  structures  and  de¬ 
tails  of  the  barotropic  tide  in  the  Korea/Tsushima  Strait. 
The  tidal  response  in  the  strait  is  strongly  frequency  de¬ 
pendent.  In  particular,  the  semi-diurnal  and  diurnal  con¬ 
stituents  exhibit  distinct  group  patterns  in  their  surface 
heights  and  currents.  Fang  and  Yang  (1988)  identified  the 
energy  inflowing  regions  of  the  strait  to  be  west  of  the 
Gotou  Islands  and  in  the  north  near  Korea  and  the 
outflowing  regions  to  be  in  the  Jeju  Strait  and  in  the  north 
near  Honshu.  In  these  various  regions  the  relative  impor¬ 
tance  of  the  semi-diurnal  and  diurnal  tides  change  with 
respect  to  each  other,  thus  representing  frequency-depend¬ 
ent  changes  induced  as  the  tides  transverse  the  strait.  Sig¬ 
nificant  differences  also  occur  in  the  basin  response  to 
Oj  compared  to  the  response  to  KP 

The  M2  tide  has  high  amplitudes  and  dominates  tidal 
surface  height  fluctuations  throughout  the  strait.  A  broad 
area  southwest  of  Tsushima  Island  exists  where  the 
semidiurnal  tide  changes  little  in  amplitude  or  phase.  The 
surface  heights  of  this  region  all  oscillate  together.  The 
model  reproduces  the  10  cm  semidiurnal  height  ampli¬ 
tude  difference  from  the  west  coast  to  the  east  coast  of 
Tsushima  Island  that  was  observed  by  Odamaki  (1989). 
Amphidromic  points  for  the  M2  and  S2  tides  appear  to  lie 
just  outside  of  the  northern  model  boundary. 

Although  the  model  semidiurnal  tide  charts  are  simi¬ 
lar  in  character  to  those  of  Odamaki  (1989),  significant 
differences  exist.  Most  notable  is  the  difference  in  am¬ 
plitude  along  the  south  coast  of  Korea,  where  our  predic¬ 
tions  are  lower  by  approximately  10  cm  than  those  of 
Odamaki  (1989).  Possibly,  the  historical  charts  are  biased 
high  because  of  their  reliance  on  coastal  data  where  tidal 
fluctuations  are  often  enhanced.  Differences  also  occur 
south  of  Jeju  Island,  but  this  region  is  near  the  model  null 
space  (Fig.  14)  and  the  model  may  not  be  able  to  predict 
the  total  tide  signal  there. 

The  diurnal  tidal  heights  are,  in  general,  smaller  than 
those  of  the  M2  and  S2  tides.  However,  in  the  northeast 
portion  of  the  strait,  the  semi-diurnal  amplitudes  rapidly 
decrease  toward  the  northeast  while  the  diurnal  ampli¬ 
tudes  remain  relatively  constant  and  thus  the  amplitudes 
of  the  Oj  and  Kx  tides  exceed  those  of  S2.  The  Kj  height 
fluctuations  are  stronger  than  those  of  Oj  for  the  half  of 
the  strait  west  of  Tsushima  Island  and  slightly  weaker  for 
the  half  east  of  Tsushima  Island.  The  model  predicts  the 
2.5-5  cm  diurnal  height  difference  from  the  east  to  the 
west  coast  of  Tsushima  Island  observed  by  Odamaki 
(1989).  Our  predictions  place  the  diurnal  amphidromic 
points  inside  the  model  domain.  Their  positions  do  not 
occur  near  the  null  space  and  therefore  should  be  accu¬ 
rate.  The  location  of  the  amphidromic  point  for  Kj  is  dif¬ 
ferent  between  our  model  and  the  charts  of  Odamaki 
(1989).  His  charts  place  the  amphidrome  against  the 


Korean  coast  just  south  of  the  city  of  Ulsan  (tidal  station 
T10),  whereas  our  model  shows  the  Kj  amphidrome  lo¬ 
cated  offshore  and  farther  southward. 

M2  currents  are  strong  throughout  the  strait  and  domi¬ 
nate  the  tidal  currents  in  most  areas.  Southwest  of 
Tsushima  Island,  the  M2  ellipses  are  weaker  and  more 
circular  than  elsewhere.  In  this  region,  the  major  semi¬ 
axes  of  the  Kx  and  the  0{  currents  are  approximately  the 
same  as  those  of  M2.  In  general,  the  strength  of  the  S2, 
K1?  and  0{  currents  are  similar.  S2  currents  are  stronger 
than  either  Kj  and  Ox  currents  north  and  south  of  Jeju 
Island,  along  the  southern  boundary,  and  near  Busan  but 
weaker  almost  everywhere  else.  For  a  broad  region  in  the 
northeast  portion  of  the  strait,  O!  and  Kj  currents  exceed 
all  other  constituents.  Here,  the  Ox  currents  are  also  sig¬ 
nificantly  larger  than  the  Kx  currents  and  thus  02  is  the 
dominate  tidal  current  constituent  for  this  region  offshore 
of  Honshu. 

Figures  11  and  12  show  the  tides  for  the  next  four 
largest  components.  Although  the  height  fluctuations  of 
N2  are  of  similar  magnitude  to  that  of  Ki  and  are  often 
greater  than  that  of  Oj,  the  tidal  currents  for  this  compo¬ 
nent  are  much  smaller  than  the  four  main  components. 
Likewise  K2  has  greater  tidal  height  fluctuations  than  Pj 
but  smaller  tidal  currents.  The  ji2  constituent  is  of  little 
significance  in  the  Korea/Tsushima  Strait. 

Kang  (1984)  and  Lee  and  Kim  (1999)  have  shown 
that  the  tides  near  Jeju  Island  exhibit  complex  local  ef¬ 
fects  that  are  difficult  to  model.  In  contrast  to  the  increas¬ 
ing  tidal  height  fluctuations  from  south  to  north  near  Jeju 
Island,  the  tidal  height  fluctuations  are  stronger  on  the 
south  side  of  the  island  than  on  the  north  side  of  the  is¬ 
land.  Thus  the  observed  M2  and  K2  amplitudes  are  70  and 
23  cm  at  station  T2,  but  are  77  and  24  cm  at  station  Tl. 
Our  model  successfully  predicts  this  increase,  giving 
amplitudes  of  69  and  23  cm  at  station  T2,  and  77  and  24 
cm  at  station  Tl.  Correctly  weighting  the  cost  function 
according  to  the  expected  error  in  the  observational 
dataset  was  required  to  obtain  this  result.  With  correct 
weighting  the  model  is  able  to  find  a  solution  that  matches 
the  observations  within  the  restrictions  of  strictly  linear 
dynamics. 

The  model  prediction  errors  (Fig.  17)  are  small  rela¬ 
tive  to  the  amplitudes  of  the  tidal  signals  in  the  Korea/ 
Tsushima  Strait  (Fig.  3)  and  are  dominated  by  bias  error 
in  the  model’s  physics,  numerics,  and/or  parameterization. 
The  model  can  only  use  linear  dynamics  to  propagate  in¬ 
coming  waves  from  the  boundary  to  match  observations. 
The  differences  shown  in  Fig.  15  indicate  the  possible 
strength  of  non-linear  effects  in  determining  the  true 
barotropic  tide  in  the  strait.  However,  it  is  likely  that  some 
portion  of  these  differences  is  not  caused  by  non-linear 
effects  but  rather  by  model  parameter  errors  such  as  er¬ 
rors  in  the  bathymetry  field.  In  particular,  bathymetric 
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sensitivity  runs  made  in  this  study  show  that  localized 
changes  in  the  model  bathymetry  can  degrade  the  model 
fit  over  broad  spatial  areas.  Odamaki  (1989)  showed  that 
the  tide-generating  force  in  the  strait  is  very  small  com¬ 
pared  to  the  other  terms  in  the  dynamical  equations.  Our 
model  agrees  with  this  finding.  If  direct  forcing  of  tides 
were  important,  differences  between  the  observations  and 
the  model  predictions  would  have  been  greater. 

Model  prediction  errors  peak  sharply  near  portions 
of  the  model  boundary  (Fig.  14).  The  model  currents  at 
these  boundaries  are  determined  by  the  observations 
through  the  inverse  method  described  in  this  paper.  How¬ 
ever,  the  error  results  do  not  suggest  that  the  forcings  so 
determined  are  invalid.  The  model  determines  with  a  high 
degree  of  accuracy  the  portion  of  the  boundary  flows  that 
impact  the  data  sites.  Additional  flows  could  exist  in  these 
regions  that  lead  to  no  impact  at  the  data  sites.  The  model 
cannot  determine  if  such  flows  exist.  Therefore,  the  high 
errors  near  the  model  boundary  represent  the  possibility 
that  the  true  barotropic  tide  could  consist  of  flows  that  do 
not  affect  the  data  sites. 

6.  Conclusions 

The  entire  Korea/Tsushima  Strait  has  high  tidal  cur¬ 
rents,  with  several  constituents  having  significant  ampli¬ 
tudes.  By  using  an  extensive  set  of  measurements  inside 
the  strait  to  predict  boundary  flows  for  a  numerical  model 
in  a  strong-constraint  inverse  problem,  tidal  predictions 
throughout  the  strait  were  obtained.  The  predictions  fol¬ 
low  strictly  linear  dynamics  and  provide  a  best  match  in 
a  least-squares  sense  to  the  observations,  weighted  ac¬ 
cording  to  their  accuracies.  Model  bias  was  estimated  by 
comparing  model-observation  differences  and  the  theo¬ 
retical  expected  values  and  ensemble  variances  of  these 
differences  calculated  only  from  the  model  dynamics  and 
measurement  accuracies.  The  total  estimated  rms  errors, 
including  model  bias,  for  M2  T]  predictions  are  less  than 
3  cm  over  most  of  the  model  domain. 

The  tides  in  the  Strait  exhibit  a  complex  structure  in 
which  diurnal  constituents  have  stronger  currents  rela¬ 
tive  to  their  sea  level  height  ranges  than  semi-diurnal 
constituents.  No  single  constituent  dominates  both  tidal 
sea  level  heights  and  velocities  in  all  parts  of  the  Strait. 
Tidal  constituents  M2,  S2,  Kj,  and  N2  have  the  high¬ 
est  ranges  in  sea  level  height.  Constituents  M2,  Kj,  S2, 
and  0{  have  the  highest  magnitudes  of  tidal  velocities. 
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